/*****************************************************************************
*
*	Calculate Mutation-Data-Matrix based on Dayhoff's procedure
*
*	Ref.	Dayhoff,M.O., Schwartz,R.M. and Orcutt,B.C. (1978)
*		"A model of evolutionary change in proteins." in Atlas of
*		Protein Sequence and Structure, vol. 5, suppl. 3, pp. 345-352
*		(Dayhoff,M.O. ed.) Natl. Biomed. Res. Found., Silver Springs.
*
*	Osamu Gotoh, ph.D.	(-2001)
*	Saitama Cancer Center Research Institute
*	818 Komuro, Ina-machi, Saitama 362-0806, Japan
*
*	Osamu Gotoh, Ph.D.	(2001-)
*	National Institute of Advanced Industrial Science and Technology
*	Computational Biology Research Center (CBRC)
*	2-41-6 Aomi, Koutou-ku, Tokyo 135-0064, Japan
*
*	Osamu Gotoh, Ph.D.      (2003-)
*	Department of Intelligence Science and Technology
*	Graduate School of Informatics, Kyoto University
*	Yoshida Honmachi, Sakyo-ku, Kyoto 606-8501, Japan
*
*	Copyright(c) Osamu Gotoh <<o.gotoh@aist.go.jp>>
*****************************************************************************/

#include <math.h>
#include "seq.h"
#include "mdm.h"

#define	DEBUG	0

#define STDSD	25.
#define GAP_WT	-60.
#define EPS	1.e-30
#define	shift_aa(x) ((x) + ALA - USE_EXG)
#define	skip_nil(x) ((x) - USE_EXG)

struct wmode {INT cmp, tab, lud, pwr;};
static int  year = 0;

static	void	usage();
#if DEBUG
static	void	printmat(FILE* fd, char* format, double matrix[][20]);
static	void	printimt(FILE* fd, char* format, double matrix[][AAS]);
#endif
static	void	ludecomp(int n, double a[][AAS]);
static	double	det(double matrix[][20]);
static	void	matcpy(double to[20][20], double from[20][20]);
static	void	matmpy(double c[20][20], double a[20][20], double b[20][20]);
static	void	matpwr(double b[20][20], double a[20][20], int n);
static	void	makes(double s[][AAS]);
static	void	putpmtx(FILE* fd, double c[20][20]);
static	void	putfmtx(FILE* fd, double c[][AAS]);
static	void	putdmtx(FILE* fd, double c[][AAS]);
static	void	pam1(double a[][20], double comp[20]);
static	void	logodd(double c[][AAS], double a[20][20], double comp[20]);
static	void	matstat(double* av, double* sd, double s[][AAS], double comp[20]);
static	double	normlz(double h[][AAS], double fact);
static	void	makmdm(FILE* fs, FILE* fd, double a[20][20], double comp[20]);
static	void	mkmdmpwr(FILE* fd, double a[20][20]);
static	void	exp_qdt(double a[20][20], double comp[20], double dt, int m);

static void usage()
{
	fputs("Usage:\n\tmakmdm [-opt...] [destpath]\n", stderr);
	fputs("Options:\n", stderr);
	fputs("\t-c-:\tSuppress output of average compositions\n", stderr);
	fputs("\t-m-:\tSuppress output of mdm matrix\n", stderr);
	fputs("\t-l:\tOutput LU-decomposited matrix\n", stderr);
	fputs("\t-p:\tOutput power matrices\n", stderr);
	fputs("\t-f78|91\n", stderr);
	fputs("\t-b:\tBrosum series\n", stderr);
	fputs("\t-v:\tVT series\n", stderr);
	fputs("If destpath is given, the files are writen there;\n", stderr);
	fputs("otherwise, in the current directry.\n", stderr);
	exit(1);
}

#if DEBUG

static void printmat(FILE* fd, char* format, double matrix[][20])
{
	for (int i = 0; i < 20; i++) {
	    for (int j = 0; j < 20; j++) {
		if (!(j % 7)) fputc('\n', fd);
		fprintf(fd, format, matrix[i][j]);
	    }
	}
	putc('\n', fd);
}

static void printimt(FILE* fd, char* format, double matrix[][AAS])
{
	for (int i = 0; i < AAS; i++) {
	    for (int j = 0; j < AAS; j++) {
		if (!(j % 12)) fputc('\n', fd);
		    fprintf(fd, format, matrix[i][j]);
	    }
	}
	putc('\n', fd);
}

#endif

/*
*	LDU decomposition of MDM
*	In results:
*
*	U[i][i] = 1		D[i][j] = 0  (i != j)
*	a[i][j] = D[i][j]  (i = j)
*	a[i][j] = U[i][j]  (i < j)
*/

static void ludecomp(int n, double a[][AAS])
{
	for (int k = 0; k < n; k++) {
	    double	p = a[k][k];
	    if (p == 0.) {
		fputs("Warning:  pivot = 0!\n", stderr);
		continue;
	    }
	    int	kp1 = k + 1;
	    for (int i = kp1; i < n; i++) {
		double	f = a[i][k] / p;
		if (f == 0.) continue;
		for (int j = kp1; j < n; j++) a[i][j] -= f * a[k][j];
	    }	/*	Normalize U */
	    for (int i = kp1; i < n; i++) a[k][i] /= p;
	}
#if DEBUG == 1
	printimt(stderr, "%8.4f ", a);
#endif
}

static double det(double matrix[][20])
{
	double	a[20][20];

	for (int j = 0; j < 20; j++)
	    for (int i = 0; i < 20; i++)
		a[j][i] = matrix[j][i];
	double dt = 1.;
	for (int k = 0; k < 20; k++) {
	    int	mx = k;
	    for (int i = k + 1; i < 20; i++)
		if (fabs(a[i][k]) > fabs(a[mx][k])) mx = i;
	    if (fabs(a[mx][k]) < EPS) return(0.);
	    if (mx != k) {
		for (int j = k; j < 20; j++) swap(a[k][j], a[mx][j]);
		if ((mx + k) % 2) dt = -dt;
	    }
	    for (int i = k + 1; i < 20; i++) {
		double	f = a[i][k]/a[k][k];
		for (int j = k + 1; j < 20; j++) a[i][j] -= f * a[k][j];
	    }
	    dt *= a[k][k];
	}
	return(dt);
}

static void matcpy(double to[20][20], double from[20][20])
{
	movmem(from, to, 400*sizeof(double));
}

static void mattr(double a[20][20])
{
	for (int i = 1; i < 20; ++i) {
	    for (int j = 0; j < i; ++j) {
		double	t = a[i][j];
		a[i][j] = a[j][i];
		a[j][i] = t;
	    }
	}
}

static void matmpy(double c[20][20], double a[20][20], double b[20][20])
{
	double	t[20][20];

	for (int j = 0; j < 20; j++) {
	    for (int i = 0; i < 20; i++) {
		t[j][i] = 0.;
		for (int k = 0; k < 20; k++)
		    t[j][i] += a[j][k]*b[k][i];
	    }
	}
	matcpy(c, t);
}

static void matpwr(double b[20][20], double a[20][20], int n)
{
	double c[20][20];

	matcpy(c, a);
	for (int i = 0; i < 20; i++) {
	    for (int j = 0; j < 20; j++) b[i][j] = 0.;
	    b[i][i] = 1.;
	}
	for ( ; n; n /= 2) {
	    if (n % 2) matmpy(b, c, b);
	    matmpy(c, c, c);
	}
}

/*
   EXTEND 20*20 MDM TO 24*24 MDM 
*/
 
static void makes(double s[][AAS])
{
	for (int i = 0; i < 20; i++) {
	    int	ii = shift_aa(i);
	    for (int j = 0; j < i; j++) {
		int	jj = shift_aa(j);
		s[jj][ii] = s[ii][jj];
	    }
	}
	for (int i = 0; i < AAS; i++) {
	    s[skip_nil(AMB)][i] = s[i][skip_nil(AMB)] = 0;
	    s[skip_nil(UNP)][i] = s[i][skip_nil(UNP)] = GAP_WT;
	}
	for (int i = 0; i < AAS; i++) {
	    s[skip_nil(ASX)][i] = s[i][skip_nil(ASX)] = 
		(s[skip_nil(ASN)][i] + s[skip_nil(ASP)][i]) / 2.;
	    s[skip_nil(GLX)][i] = s[i][skip_nil(GLX)] = 
		(s[skip_nil(GLN)][i] + s[skip_nil(GLU)][i]) / 2.;
	}
 	s[skip_nil(UNP)][skip_nil(UNP)] = 0;
	s[skip_nil(AMB)][skip_nil(AMB)] = 1.;
}

static void putpmtx(FILE* fd, double c[20][20])
{
	fwrite(c, sizeof(double), 20 * 20, fd);
}

static void putfmtx(FILE* fd, double c[][AAS])
{
	double	s[AASCMB];

	int	k = 0;
	for (int i = 0; i < AAS; i++)
	    for (int j = 0; j <= i; ) 
		s[k++] = c[i][j++];
	fwrite(s, sizeof(double), AASCMB, fd);
}

static void putdmtx(FILE* fd, double c[][AAS])
{
	double	b[AASCMB];

	ludecomp(AAS, c);
	int	k = 0;
	for (int i = 0; i < AAS; i++)
	    for (int j = 0; j <= i; )
		b[k++] = c[j++][i];
	fwrite(b, sizeof(double), AASCMB, fd);
}

/*   MUTATION PROBABILITY MATRIX AT 1 PAM	*/

static void pam1(double a[][20], double comp[20])
{ 
static	double rmt[2][20] = {
	    {100.,83.,104.,86.,44.,84.,77.,50.,91.,103.,54.,
		72.,93.,51.,58.,117.,107.,25.,50.,98.}, 
	    {100.,65.,134.,106.,20.,93.,102.,49.,66.,96.,40.,
		56.,94.,41.,56.,120.,97.,18.,41.,74.}
	};
static	int rawdat[2][190] = {
	{
	247,
	216,116,
	386,48,1433,
	106,125,32,13,
	208,750,159,130,9,
	600,119,180,2914,8,1027,
	1183,614,291,577,98,84,610,
	46,446,466,144,40,635,41,41,
	173,76,130,37,19,20,43,25,26,
	257,205,63,34,36,314,65,56,134,1324,
	200,2348,758,102,7,858,754,142,85,75,94,
	100,61,39,27,23,52,30,27,21,704,974,103,
	51,16,15,8,66,9,13,18,50,196,1093,7,49,
	901,217,31,39,15,395,71,93,157,31,578,77,23,36,
	2413,413,1738,244,353,182,156,1131,138,172,436,228,54,309,1138,
	2440,230,693,151,66,149,142,164,76,930,172,398,343,39,412,2258,
	11,109,2,5,38,12,12,69,5,12,82,9,8,37,6,36,8,
	41,46,114,89,164,40,15,15,514,61,84,20,17,850,22,164,45,41,
	1766,69,55,127,99,58,226,276,22,3938,1261,58,559,189,84,219,526,27,42
	}, {
	30,
	109,17,
	154,0,532,
	33,10,0,0,
	93,120,50,76,0,
	266,0,94,831,0,422,
	579,10,156,162,10,30,112,
	21,103,226,43,10,243,23,10,
	66,30,36,13,17,8,35,0,3,
	95,17,37,0,0,75,15,17,40,253,
	57,477,322,85,0,147,104,60,23,43,39,
	29,17,0,0,0,20,7,7,0,57,207,90,
	20,7,7,0,0,0,0,17,20,90,167,0,17,
	345,67,27,10,10,93,40,49,50,7,43,43,4,7,
	772,137,432,98,117,47,86,450,26,20,32,168,20,40,269,
	590,20,169,57,10,37,31,50,14,129,52,200,28,10,73,696,
	0,27,3,0,0,0,0,0,3,0,13,0,0,10,0,17,0,
	20,3,36,0,30,0,10,0,40,13,23,10,0,260,0,22,23,6,
	365,20,13,17,33,27,37,97,30,661,303,17,77,10,50,43,186,0,17
	}
	};
	double b[20][20];
	double count[20][20];
	double	delta = 0.01;

	int	k = 0;
	for (int i = 0; i < 20; i++) {
	    for (int j = 0; j < i; j++)
		count[i][j] = count[j][i] = rawdat[year][k++];
	    count[i][i] = 0.;
	}
	for (int j = 0; j < 20; j++) {
	    double	sum = 0;
	    for (int i = 0; i < 20; i++)
		sum += count[i][j];
	    for (int i = 0; i < 20; i++)
		a[i][j] = delta * rmt[year][j] * count[i][j] / sum;
	    a[j][j] = -delta * rmt[year][j];
	}
	matcpy(b, a);

	double	dt = 0;
	double	sum = 0;
	for (int i = 0; i < 20; i++) {
	    for (int j = 0; j < 20; j++) {
		if (i > 0) b[i-1][j] = a[i-1][j];
		b[i][j] = 0.;
	    }
	    b[i][i] = 1.;
	    comp[i] = det(b);
	    dt += comp[i];
	    sum += comp[i] * a[i][i];
	}

	sum = -0.01 * dt / sum;
	for (int i = 0; i < 20; i++) {
	    comp[i] /= dt;
	    for (int j = 0; j < 20; j++) a[i][j] *= sum;
	    a[i][i] += 1.;
	}
#if DEBUG == 2
	printmat(stdout, "%8.5f ", a);
#endif
}


/*   MUTATION DATA MATRIX AT NN PAM	*/


static void logodd(double c[][AAS], double a[20][20], double comp[20])
{
	for (int i = 0; i < 20; i++) {
	    int	ii = shift_aa(i);
	    for (int j = 0; j <= i; j++) 
		c[ii][shift_aa(j)] = log(a[i][j] / comp[i]);
	}
}

/*   MEAN VALUE OF 'S' FOR A GIVEN PAIR OF COMPOSITIONS	*/

static void matstat(double* av, double* sd, double s[][AAS], double comp[20])
{
	*av = *sd = 0.;
	for (int i = 0; i < 20; i++) {
	    double	aw = 0.;
	    double	sw = 0.;
	    int	ii = shift_aa(i);
	    for (int j = 0; j < 20; j++) {
		int	jj = shift_aa(j);
		double	t = i > j? s[ii][jj]: s[jj][ii];
		aw += t * comp[j];
		sw += t * t * comp[j];
	    }
	    *av += aw * comp[i];
	    *sd += sw * comp[i];
	}
	*sd = sqrt(*sd - *av * *av);
}

static void exp_qdt(double a[20][20], double comp[20], double dt, int m)
{
static double q[2][20][20] = {
{
  { -0.01141100,     0.00031700,     0.00029000,     0.00037600,     0.00018200, 
     0.00046200,     0.00075300,     0.00111500,     0.00010700,     0.00027700, 
     0.00054000,     0.00045000,     0.00021300,     0.00014600,     0.00066600, 
     0.00230100,     0.00152000,     0.00003600,     0.00009100,     0.00156900},
  {  0.00048500,    -0.00943100,     0.00037200,     0.00020900,     0.00010700, 
     0.00114600,     0.00043500,     0.00043800,     0.00054600,     0.00013400, 
     0.00041700,     0.00300500,     0.00013900,     0.00008900,     0.00031800, 
     0.00069200,     0.00038600,     0.00009700,     0.00014800,     0.00026800},
  {  0.00047100,     0.00039800,    -0.01172600,     0.00214800,     0.00007100, 
     0.00062400,     0.00064300,     0.00080300,     0.00069500,     0.00027600, 
     0.00020200,     0.00133200,     0.00010300,     0.00010500,     0.00016800, 
     0.00221700,     0.00112300,     0.00001200,     0.00018000,     0.00015500},
  {  0.00054000,     0.00019000,     0.00184800,    -0.00948400,     0.00001600, 
     0.00034100,     0.00312800,     0.00064100,     0.00024000,     0.00006600, 
     0.00017700,     0.00043600,     0.00004400,     0.00003400,     0.00025700, 
     0.00079600,     0.00043500,     0.00001600,     0.00012200,     0.00015700},
  {  0.00096900,     0.00036300,     0.00023300,     0.00006400,    -0.00712500, 
     0.00009700,     0.00015300,     0.00036700,     0.00016800,     0.00028700, 
     0.00070800,     0.00010400,     0.00013700,     0.00035300,     0.00010500, 
     0.00129100,     0.00050900,     0.00010400,     0.00034700,     0.00076600},
  {  0.00086600,     0.00140400,     0.00071200,     0.00045700,     0.00003300, 
    -0.01220200,     0.00218000,     0.00029300,     0.00088600,     0.00015300, 
     0.00066000,     0.00182200,     0.00021300,     0.00007500,     0.00056000, 
     0.00083500,     0.00062400,     0.00005800,     0.00015300,     0.00021800},
  {  0.00091700,     0.00034600,     0.00047500,     0.00264400,     0.00003100, 
     0.00142600,    -0.01022200,     0.00046200,     0.00018200,     0.00012300, 
     0.00027800,     0.00130300,     0.00008300,     0.00004700,     0.00036300, 
     0.00065800,     0.00045900,     0.00002700,     0.00009500,     0.00030300},
  {  0.00130600,     0.00033400,     0.00056800,     0.00052100,     0.00008100, 
     0.00018200,     0.00044700,    -0.00603900,     0.00009500,     0.00009000, 
     0.00018000,     0.00027100,     0.00006600,     0.00007700,     0.00018700, 
     0.00110000,     0.00020800,     0.00006000,     0.00005600,     0.00021000},
  {  0.00037500,     0.00125800,     0.00144600,     0.00058500,     0.00011300, 
     0.00166900,     0.00052700,     0.00028500,    -0.01107500,     0.00017800, 
     0.00053500,     0.00054300,     0.00015300,     0.00028800,     0.00044800, 
     0.00070800,     0.00046800,     0.00004100,     0.00127200,     0.00018300},
  {  0.00033700,     0.00010700,     0.00021600,     0.00005800,     0.00006400, 
     0.00009700,     0.00012800,     0.00009300,     0.00006300,    -0.01193400, 
     0.00297500,     0.00018600,     0.00076700,     0.00047900,     0.00011300, 
     0.00021600,     0.00076800,     0.00003500,     0.00014000,     0.00509200},
  {  0.00041800,     0.00021300,     0.00009700,     0.00009800,     0.00010200, 
     0.00027700,     0.00018100,     0.00011800,     0.00012000,     0.00181000, 
    -0.00779500,     0.00016800,     0.00090300,     0.00091800,     0.00024900, 
     0.00032600,     0.00029300,     0.00009800,     0.00020400,     0.00120200},
  {  0.00058400,     0.00255900,     0.00103000,     0.00039900,     0.00002400, 
     0.00125800,     0.00139500,     0.00029700,     0.00020200,     0.00019000, 
     0.00028400,    -0.01070600,     0.00019200,     0.00008200,     0.00033400, 
     0.00071100,     0.00073300,     0.00003800,     0.00013000,     0.00026400},
  {  0.00074700,     0.00032100,     0.00022500,     0.00011500,     0.00009000, 
     0.00039900,     0.00024400,     0.00020100,     0.00015500,     0.00201200, 
     0.00397900,     0.00051600,    -0.01262900,     0.00049600,     0.00012000, 
     0.00038800,     0.00096900,     0.00009700,     0.00018200,     0.00137300},
  {  0.00026700,     0.00010700,     0.00012000,     0.00004700,     0.00012600, 
     0.00007400,     0.00007700,     0.00012400,     0.00015100,     0.00069100, 
     0.00216100,     0.00011800,     0.00025600,    -0.00751400,     0.00009000, 
     0.00040300,     0.00019900,     0.00021500,     0.00170800,     0.00058000},
  {  0.00106700,     0.00033200,     0.00017100,     0.00029700,     0.00003200, 
     0.00047700,     0.00049300,     0.00026000,     0.00020300,     0.00015000, 
     0.00049200,     0.00042300,     0.00005500,     0.00007900,    -0.00667200, 
     0.00126600,     0.00056200,     0.00002800,     0.00008000,     0.00020500},
  {  0.00252500,     0.00048800,     0.00143000,     0.00060700,     0.00027100, 
     0.00048100,     0.00058800,     0.00102300,     0.00021900,     0.00018300, 
     0.00044000,     0.00059600,     0.00011800,     0.00023000,     0.00086600, 
    -0.01313500,     0.00255400,     0.00004300,     0.00015700,     0.00031600},
  {  0.00208200,     0.00034000,     0.00091200,     0.00041200,     0.00012800, 
     0.00044700,     0.00050900,     0.00023700,     0.00017900,     0.00078500, 
     0.00050800,     0.00077200,     0.00037800,     0.00014400,     0.00047300, 
     0.00318200,    -0.01282300,     0.00002900,     0.00014200,     0.00116400},
  {  0.00021300,     0.00038900,     0.00004200,     0.00006800,     0.00012100, 
     0.00018500,     0.00013500,     0.00031200,     0.00006900,     0.00017200, 
     0.00075300,     0.00017500,     0.00016200,     0.00070500,     0.00010300, 
     0.00024600,     0.00012900,    -0.00491800,     0.00074800,     0.00019100},
  {  0.00020800,     0.00022500,     0.00025500,     0.00020500,     0.00016200, 
     0.00018800,     0.00019000,     0.00011200,     0.00086700,     0.00026400, 
     0.00062000,     0.00023400,     0.00012100,     0.00218500,     0.00011500, 
     0.00034500,     0.00024600,     0.00029200,    -0.00717300,     0.00033900},
  {  0.00183700,     0.00019700,     0.00010900,     0.00012500,     0.00016700, 
     0.00013200,     0.00028900,     0.00020700,     0.00005900,     0.00444900, 
     0.00176700,     0.00023600,     0.00045800,     0.00035800,     0.00014400, 
     0.00033600,     0.00099000,     0.00003700,     0.00016500,    -0.01206200},
 }, {
  { -0.01098000,     0.00044400,     0.00021600,     0.00021800,     0.00024900, 
     0.00038000,     0.00054400,     0.00129200,     0.00016400,     0.00037500, 
     0.00076600,     0.00050500,     0.00038000,     0.00023600,     0.00050200, 
     0.00219600,     0.00075900,     0.00004500,     0.00021000,     0.00149900},
  {  0.00062400,    -0.01044200,     0.00053000,     0.00031100,     0.00008400, 
     0.00105300,     0.00083100,     0.00045200,     0.00035700,     0.00025000, 
     0.00055900,     0.00267500,     0.00020100,     0.00013500,     0.00027400, 
     0.00073300,     0.00063700,     0.00011800,     0.00026400,     0.00035400},
  {  0.00043400,     0.00075700,    -0.01156700,     0.00146800,     0.00012400, 
     0.00072600,     0.00073500,     0.00109200,     0.00065600,     0.00026600, 
     0.00033600,     0.00096100,     0.00017800,     0.00021100,     0.00034200, 
     0.00157000,     0.00094600,     0.00005100,     0.00029600,     0.00041800},
  {  0.00036800,     0.00037400,     0.00123700,    -0.00888500,     0.00008700, 
     0.00054500,     0.00191400,     0.00073700,     0.00023000,     0.00013000, 
     0.00027600,     0.00057400,     0.00004300,     0.00013200,     0.00031000, 
     0.00099200,     0.00049600,     0.00004600,     0.00016600,     0.00022800},
  {  0.00066000,     0.00015800,     0.00016400,     0.00013600,    -0.00492900, 
     0.00011000,     0.00009200,     0.00028400,     0.00015000,     0.00029700, 
     0.00047600,     0.00016800,     0.00013200,     0.00030700,     0.00011200, 
     0.00054500,     0.00036600,     0.00004800,     0.00018400,     0.00054000},
  {  0.00084600,     0.00167000,     0.00080700,     0.00071800,     0.00009200, 
    -0.01442900,     0.00225000,     0.00067300,     0.00050800,     0.00021100, 
     0.00073900,     0.00178200,     0.00051500,     0.00018300,     0.00039100, 
     0.00115500,     0.00091400,     0.00013500,     0.00033500,     0.00050500},
  {  0.00076900,     0.00083500,     0.00051800,     0.00160000,     0.00004900, 
     0.00142600,    -0.01073000,     0.00043300,     0.00031200,     0.00023700, 
     0.00039800,     0.00122200,     0.00013000,     0.00011700,     0.00052500, 
     0.00079900,     0.00065000,     0.00007200,     0.00019000,     0.00044800},
  {  0.00125100,     0.00031100,     0.00052700,     0.00042200,     0.00010400, 
     0.00029300,     0.00029600,    -0.00618800,     0.00016100,     0.00010900, 
     0.00029100,     0.00032500,     0.00008900,     0.00018900,     0.00026100, 
     0.00088600,     0.00029100,     0.00005500,     0.00012600,     0.00020100},
  {  0.00041200,     0.00063900,     0.00082300,     0.00034300,     0.00014200, 
     0.00057300,     0.00055500,     0.00041800,    -0.00779300,     0.00016400, 
     0.00033700,     0.00055100,     0.00011200,     0.00024000,     0.00018400, 
     0.00061400,     0.00061800,     0.00009300,     0.00070300,     0.00027200},
  {  0.00047300,     0.00022400,     0.00016700,     0.00009700,     0.00014100, 
     0.00012000,     0.00021200,     0.00014200,     0.00008200,    -0.01327600, 
     0.00343200,     0.00025800,     0.00080400,     0.00060900,     0.00014300, 
     0.00028800,     0.00059000,     0.00006900,     0.00023900,     0.00518600},
  {  0.00060400,     0.00031400,     0.00013200,     0.00012900,     0.00014200, 
     0.00026100,     0.00022200,     0.00023700,     0.00010600,     0.00214500, 
    -0.00948000,     0.00025700,     0.00120600,     0.00101500,     0.00018400, 
     0.00024200,     0.00041800,     0.00019200,     0.00027100,     0.00140300},
  {  0.00073400,     0.00276600,     0.00069700,     0.00049300,     0.00009200, 
     0.00116200,     0.00125700,     0.00048800,     0.00031900,     0.00029800, 
     0.00047300,    -0.01182100,     0.00021800,     0.00013300,     0.00042400, 
     0.00084600,     0.00072400,     0.00007700,     0.00023200,     0.00038800},
  {  0.00130900,     0.00049400,     0.00030600,     0.00008900,     0.00017200, 
     0.00079700,     0.00031700,     0.00031700,     0.00015400,     0.00219700, 
     0.00527100,     0.00051800,    -0.01652400,     0.00094100,     0.00019500, 
     0.00059300,     0.00085200,     0.00012700,     0.00036600,     0.00150900},
  {  0.00039600,     0.00016100,     0.00017700,     0.00013100,     0.00019400, 
     0.00013800,     0.00013900,     0.00032700,     0.00016000,     0.00081000, 
     0.00216000,     0.00015300,     0.00045800,    -0.00926300,     0.00021400, 
     0.00044700,     0.00032700,     0.00035400,     0.00185900,     0.00065800},
  {  0.00090300,     0.00035000,     0.00030700,     0.00032900,     0.00007600, 
     0.00031600,     0.00066700,     0.00048500,     0.00013100,     0.00020300, 
     0.00041900,     0.00052500,     0.00010100,     0.00022900,    -0.00733900, 
     0.00094000,     0.00058900,     0.00008700,     0.00019300,     0.00048900},
  {  0.00243200,     0.00057800,     0.00086800,     0.00065000,     0.00022800, 
     0.00057400,     0.00062600,     0.00101300,     0.00027100,     0.00025300, 
     0.00034000,     0.00064500,     0.00019100,     0.00029500,     0.00057900, 
    -0.01310700,     0.00275000,     0.00009900,     0.00025300,     0.00046200},
  {  0.00101600,     0.00060800,     0.00063200,     0.00039300,     0.00018500, 
     0.00054900,     0.00061600,     0.00040200,     0.00032900,     0.00062600, 
     0.00071000,     0.00066700,     0.00033100,     0.00026100,     0.00043900, 
     0.00332500,    -0.01288000,     0.00013700,     0.00023800,     0.00141600},
  {  0.00021500,     0.00040400,     0.00012300,     0.00013100,     0.00008700, 
     0.00029200,     0.00024300,     0.00027200,     0.00017700,     0.00026300, 
     0.00116700,     0.00025600,     0.00017700,     0.00101200,     0.00023300, 
     0.00043000,     0.00049000,    -0.00779800,     0.00138100,     0.00044500},
  {  0.00044100,     0.00039400,     0.00031100,     0.00020700,     0.00014600, 
     0.00031600,     0.00028200,     0.00027400,     0.00058700,     0.00039800, 
     0.00072200,     0.00033600,     0.00022300,     0.00232600,     0.00022600, 
     0.00047900,     0.00037300,     0.00060400,    -0.00908000,     0.00043500},
  {  0.00158500,     0.00026700,     0.00022100,     0.00014300,     0.00021500, 
     0.00024000,     0.00033500,     0.00022000,     0.00011400,     0.00434500, 
     0.00188100,     0.00028200,     0.00046300,     0.00041500,     0.00028700, 
     0.00044100,     0.00111800,     0.00009800,     0.00021900,    -0.01288900},
 }
};

static double u[2][20][20] = {
 {
  {  0.25587860,     0.01247304,     0.01614590,     0.00129445,     0.00813324, 
    -0.00375616,     0.00299541,    -0.00193491,    -0.00046507,     0.62925300, 
    -0.06283629,    -0.00658223,    -0.03831081,     0.00114387,     0.00064285, 
    -0.19847435,     0.08668531,     0.00063772,     0.00358894,    -0.70651253},
  { -0.13490481,    -0.03637193,    -0.17027114,     0.03220705,    -0.01084565, 
     0.00105309,    -0.02173227,    -0.03991957,     0.01268758,     0.21721438, 
    -0.05723474,     0.06195394,     0.02394371,    -0.01021990,    -0.03497485, 
     0.79671424,    -0.51456058,    -0.00044032,    -0.00086110,    -0.11343715},
  { -0.32660494,    -0.15583175,    -0.50665984,     0.38046780,     0.00392883, 
     0.25957352,    -0.41168724,     0.05692845,     0.02902241,     0.03478155, 
     0.02985005,     0.22151636,    -0.09440714,     0.00205976,    -0.01462007, 
     0.07242858,     0.45481475,     0.00009772,    -0.00562122,    -0.03003758},
  { -0.49578789,     0.05972667,    -0.14727153,    -0.20874113,    -0.00351799, 
    -0.53853350,     0.46948433,     0.04435550,     0.15462261,     0.06667100, 
     0.05345465,     0.05327601,    -0.06905852,    -0.00410851,     0.00442621, 
     0.21785861,     0.43307071,     0.00138949,    -0.02430252,    -0.06701419},
  {  0.51920670,    -0.26908058,    -0.16708941,     0.14325383,     0.00226966, 
    -0.41774062,    -0.11942581,    -0.02708011,     0.21062598,    -0.12037574, 
     0.03114495,     0.61602979,     0.05224574,     0.00564441,     0.00257203, 
    -0.22306825,    -0.24083005,     0.00269229,    -0.03268529,     0.03169046},
  { -0.15837064,     0.06624980,    -0.16361502,     0.10368773,     0.00481378, 
    -0.03260845,    -0.01177384,     0.02604392,     0.02089793,    -0.05128667, 
    -0.55643399,    -0.12315282,     0.87860642,     0.03575674,     0.02403161, 
    -0.05310337,     0.17276931,    -0.00317630,    -0.01212541,    -0.16721072},
  {  0.42608788,     0.51627082,    -0.35590198,     0.31766880,    -0.01758642, 
    -0.18104552,    -0.05756374,    -0.06026914,     0.12667808,    -0.10967947, 
     0.08964577,    -0.55987438,    -0.06505238,     0.00203453,    -0.03426199, 
     0.01436271,    -0.04415511,    -0.00318037,    -0.02356448,     0.01938639},
  {  0.24639074,    -0.30270098,    -0.35431646,    -0.42325891,    -0.00011715, 
     0.49297724,     0.49537675,     0.02312025,     0.37622598,    -0.04291207, 
    -0.09814598,    -0.18058854,    -0.00821367,     0.07239527,    -0.06929580, 
    -0.09254076,    -0.04970006,     0.00265090,    -0.12073049,     0.03338374},
  {  0.46211259,    -0.04270917,    -0.42481336,    -0.18576411,    -0.01385986, 
     0.02863488,     0.43061922,    -0.05576888,    -0.47887883,    -0.31564803, 
     0.37591398,     0.17185859,     0.04866232,    -0.14611551,    -0.09078102, 
     0.14243205,     0.13869476,    -0.00944315,     0.21003509,    -0.24518153},
  {  0.09899555,    -0.13565822,     0.27859624,    -0.07944823,    -0.02694513, 
     0.02481813,    -0.29811602,    -0.09114343,     0.18980855,    -0.36749757, 
     0.63712537,    -0.12816231,     0.07510161,    -0.12453481,    -0.10803033, 
     0.32876446,     0.20687504,    -0.00049682,    -0.02824895,    -0.45180314},
  { -0.40174026,     0.02069862,    -0.09219145,     0.13845354,     0.05886104, 
     0.05252039,     0.14320370,     0.19441338,     0.11847455,     0.07946834, 
     0.47890919,    -0.04619009,     0.05913996,    -0.54385752,     0.15217489, 
    -0.36203509,    -0.33151190,    -0.00568782,     0.31938655,    -0.03249000},
  { -0.35996508,     0.11452579,    -0.07082616,     0.17389034,     0.08826188, 
     0.04776611,     0.20434056,     0.15333291,    -0.07525246,    -0.21322951, 
     0.51365507,     0.08653449,     0.01880013,     0.43599755,     0.14652420, 
    -0.24545392,    -0.30250048,    -0.00745756,    -0.36705682,    -0.34188702},
  {  0.04977645,    -0.59641256,     0.12778944,     0.56003424,     0.03407241, 
    -0.09969812,     0.36248657,    -0.19039296,    -0.05574943,     0.02294632, 
    -0.02929363,    -0.38391603,    -0.01552149,     0.03310886,     0.02227762, 
     0.05587035,     0.04138841,     0.01814335,     0.00793253,     0.03515768},
  { -0.07425982,    -0.12189637,    -0.02001455,     0.08289260,    -0.90825813, 
    -0.00760918,     0.05867650,     0.37257146,    -0.01838309,     0.08563510, 
     0.26033957,    -0.06391674,     0.03335200,     0.15784684,     0.30361986, 
    -0.11694552,    -0.08165021,     0.00382861,     0.03989572,     0.01427536},
  { -0.11158234,     0.16842033,     0.10730048,     0.23781995,    -0.09476092, 
     0.08676714,     0.23449922,    -0.05303636,     0.02069676,     0.10937553, 
     0.13703575,     0.19720411,     0.04019537,    -0.07311724,    -0.87585337, 
    -0.10573429,    -0.02935961,     0.01255802,    -0.09422493,     0.08579640},
  { -0.06074292,     0.19410111,     0.01778206,     0.06760008,    -0.07270542, 
     0.12756888,     0.15297036,    -0.90714117,     0.02644123,     0.06254747, 
     0.09123925,     0.21647182,     0.02560326,    -0.17686561,     0.27630412, 
    -0.02578959,     0.04565278,     0.06058316,    -0.17010592,     0.04848506},
  {  0.18547068,    -0.14840423,    -0.08439249,    -0.14405112,     0.01673944, 
    -0.09742666,    -0.14306890,     0.34859566,    -0.14764417,     0.31869711, 
     0.34199894,    -0.15283393,     0.08113362,    -0.44814494,     0.05755438, 
     0.04761285,     0.07891294,     0.18369580,    -0.62971105,     0.33526607},
  {  0.08212606,    -0.13574213,    -0.07574361,    -0.14267950,     0.01464079, 
    -0.06671490,    -0.13270069,    -0.24380218,    -0.04980758,     0.47022341, 
     0.67011039,    -0.11689276,     0.12750416,     0.03945840,    -0.00301924, 
    -0.02025950,     0.07769739,    -0.78107460,    -0.16390177,     0.45057785},
  {  0.12701046,     0.17286080,     0.18525368,     0.29966933,    -0.02859237, 
     0.13585434,     0.30423494,     0.36740223,     0.03020688,    -0.33542331, 
    -0.64148038,     0.22790784,    -0.10207524,    -0.41452435,     0.17019225, 
     0.18158819,     0.07178403,    -0.20612860,    -0.26409692,    -0.28164373},
  {  0.34468347,     0.22396080,     0.20681205,     0.24066721,     0.06459242, 
     0.18288345,     0.28338247,     0.29335133,     0.09786315,     0.26472097, 
     0.43661849,     0.26479548,     0.09869082,     0.18479316,     0.21333349, 
     0.31638782,     0.25392472,     0.05667884,     0.14471526,     0.29928055}
 }, {
  { -0.16780967,     0.00148279,    -0.00893060,    -0.00162317,    -0.00363427, 
    -0.05215719,     0.00748901,     0.00474903,     0.00100110,    -0.70823713, 
     0.03862533,     0.00948079,     0.18682470,    -0.00312308,    -0.00570827, 
     0.13305495,    -0.14076313,    -0.00295770,    -0.00042293,     0.71265941},
  {  0.01904986,    -0.01158068,     0.00445472,    -0.01091234,    -0.00052188, 
     0.17599220,    -0.04807103,    -0.00417274,    -0.00386593,    -0.27666307, 
     0.56005664,     0.00152358,    -0.83581097,     0.02245814,    -0.00752789, 
     0.01153128,     0.00673393,    -0.00552626,     0.00586346,     0.39698896},
  { -0.22700082,     0.00197043,    -0.05292338,    -0.02320830,    -0.00121402, 
     0.00365632,     0.00904964,    -0.01947804,     0.00557599,     0.21001074, 
    -0.00945560,     0.00544237,    -0.03552129,    -0.01016888,    -0.00764081, 
     0.77247750,    -0.56002175,     0.00183916,    -0.00054879,    -0.06284048},
  {  0.02750815,     0.09465889,     0.07214392,    -0.03485309,    -0.00025952, 
    -0.90109686,     0.33544047,     0.01373128,     0.02059211,     0.00578438, 
     0.15132139,     0.22039989,    -0.11261137,    -0.00938646,    -0.00573528, 
     0.06257254,     0.07555040,     0.00373930,     0.01743728,    -0.03693741},
  { -0.01546003,    -0.56973920,    -0.09040811,     0.03759865,    -0.00048763, 
     0.03669543,    -0.17990119,     0.00411987,     0.00710832,     0.00402534, 
    -0.00536652,     0.79582733,     0.00952610,     0.00001084,    -0.00979796, 
    -0.01810342,    -0.00018181,     0.00303727,     0.00151559,    -0.01001889},
  { -0.61332283,     0.07243145,    -0.47946882,     0.10462059,     0.00221052, 
    -0.01082619,    -0.03811263,     0.10425331,     0.01898426,     0.01292363, 
     0.04697540,    -0.04181634,    -0.00483410,    -0.02597590,     0.00111657, 
     0.29643217,     0.59167724,    -0.01541346,     0.01961159,    -0.04146647},
  { -0.60457850,    -0.28364424,     0.46655902,    -0.38358347,     0.00772344, 
     0.10934629,     0.38931405,     0.05570607,    -0.06290997,     0.11377160, 
    -0.13447261,    -0.03125621,    -0.01290422,     0.03289816,    -0.00908017, 
     0.04991002,     0.23335316,    -0.00148675,    -0.02258756,     0.08792189},
  { -0.07261727,     0.11647976,    -0.08663435,     0.07745023,    -0.00210836, 
    -0.04880187,    -0.03251894,     0.04218721,     0.03222895,     0.32417246, 
    -0.77391109,     0.05769854,    -0.09267894,     0.32100450,     0.01189640, 
    -0.15037376,    -0.05463458,     0.02188663,    -0.17723066,     0.48650517},
  {  0.31148493,    -0.22873208,    -0.44729578,    -0.18687763,     0.00289954, 
     0.12279095,     0.76596681,     0.01717810,     0.04835641,    -0.07247408, 
    -0.01119278,    -0.05735444,     0.00365662,     0.09562966,    -0.04574027, 
    -0.09190875,    -0.13677628,     0.00362798,    -0.05642514,    -0.03681377},
  {  0.08656048,    -0.00568199,     0.09440300,    -0.03054189,    -0.00847183, 
    -0.00240766,    -0.10612810,    -0.08386034,     0.04079078,    -0.29603620, 
     0.23891930,    -0.00728584,     0.02878084,     0.60989390,    -0.02394822, 
     0.23018567,     0.13434821,     0.02621993,    -0.55555679,    -0.37018320},
  {  0.54642641,    -0.25313359,     0.04956620,    -0.09455112,    -0.03629426, 
    -0.03809157,    -0.10064688,    -0.23309880,    -0.05506002,    -0.05325204, 
    -0.41978463,    -0.16629162,    -0.03539840,    -0.03657718,    -0.14239382, 
     0.48090840,     0.36660431,    -0.00851244,     0.12702947,     0.10255160},
  {  0.09425674,     0.59686608,    -0.09337933,    -0.61128012,    -0.00190290, 
     0.09082771,    -0.08512037,     0.00666671,    -0.18841809,    -0.12648110, 
    -0.14989550,     0.41104069,    -0.00473447,     0.12485610,    -0.06438639, 
     0.04100955,     0.03056000,    -0.12398935,     0.17286632,    -0.11936233},
  { -0.15186556,    -0.18017376,     0.07179163,     0.47376969,     0.00222149, 
    -0.01626849,     0.14917222,    -0.02426529,    -0.29808725,    -0.13499240, 
    -0.06957519,    -0.10031750,    -0.01680502,     0.59532110,     0.01613924, 
    -0.05673906,    -0.10485010,    -0.44001416,     0.47854983,    -0.19301141},
  { -0.02351994,     0.18895716,     0.03813591,     0.38975699,     0.02117132, 
     0.07777915,     0.26105050,    -0.10559074,    -0.76160131,    -0.03649882, 
    -0.02521953,     0.17057628,    -0.00343131,    -0.08400831,    -0.35676787, 
     0.03891838,     0.02501690,     0.32407513,    -0.10381314,    -0.03498673},
  { -0.04038258,     0.14221912,     0.08851927,     0.13551665,    -0.01370165, 
     0.05719616,     0.07885548,     0.02012803,     0.28121993,     0.12005879, 
     0.15605872,     0.10092953,     0.03174647,    -0.07691265,    -0.92090244, 
    -0.02474205,     0.01650050,    -0.16972067,    -0.08181924,     0.09923264},
  {  0.19503329,     0.02526760,    -0.02747247,     0.01243242,    -0.02408995, 
     0.01423852,     0.06940205,    -0.30940010,    -0.31283032,     0.36089567, 
     0.37353465,     0.04860153,     0.07118211,    -0.37991630,     0.25365023, 
     0.07773656,     0.11326489,    -0.42092651,    -0.54391815,     0.40331424},
  {  0.12080271,    -0.21187801,    -0.07508127,    -0.18278705,    -0.03108173, 
    -0.09288451,    -0.21618692,     0.89607889,    -0.19809083,     0.10532456, 
     0.17281740,    -0.18949249,     0.02390394,     0.05460723,    -0.13431498, 
    -0.02432170,    -0.05385047,    -0.01951479,    -0.04298244,     0.09893248},
  { -0.06938998,    -0.11482077,    -0.11907196,    -0.21164564,    -0.26642050, 
    -0.06343515,    -0.15710491,    -0.53334001,    -0.06774907,     0.35980187, 
     0.60787674,    -0.11996318,     0.09688977,     0.30914090,    -0.10078750, 
    -0.13629889,    -0.02544865,     0.11179607,     0.18120963,     0.31876127},
  {  0.04360046,     0.13957694,     0.09792894,     0.16892399,    -0.88683678, 
     0.07753688,     0.16642347,     0.29624438,     0.06576026,    -0.11139253, 
    -0.17376742,     0.13454130,    -0.02429519,    -0.09694898,     0.11562334, 
     0.09135183,     0.04699879,    -0.01408698,    -0.03338638,    -0.10379632},
  { -0.33781904,    -0.24044139,    -0.16853223,    -0.19992065,    -0.12747441, 
    -0.15164179,    -0.23907228,    -0.34896733,    -0.13436519,    -0.26775500, 
    -0.42835608,    -0.23249305,    -0.09798896,    -0.20133370,    -0.18798585, 
    -0.30499824,    -0.25226955,    -0.07037249,    -0.16083718,    -0.31951152},
 }
};

static double v[2][20][20] = {
 {
  {  0.20231005,    -0.11211708,    -0.20926420,    -0.30846361,     0.31981624, 
    -0.05136874,     0.27737239,     0.12191684,     0.20937843,     0.06865873, 
    -0.23673655,    -0.21376173,     0.03165520,    -0.01540672,    -0.06783700, 
    -0.04242846,     0.07940551,     0.01641931,     0.06138131,     0.22360680},
  {  0.01602655,    -0.04583743,    -0.16004046,     0.06112554,    -0.25454621, 
     0.03406625,     0.51491364,    -0.22432170,    -0.03111967,    -0.14516803, 
     0.01778461,     0.10382391,    -0.61946639,    -0.04050668,     0.15136117, 
     0.20306324,    -0.09903340,    -0.04154125,     0.12759882,     0.22360680},
  {  0.02260681,    -0.24095280,    -0.56854921,    -0.15725676,    -0.18356471, 
    -0.08992730,    -0.39755405,    -0.29426521,    -0.32954636,     0.32588005, 
    -0.09369993,    -0.07249983,     0.14531659,    -0.00748082,     0.10712145, 
     0.02092995,    -0.06080510,    -0.02566303,     0.14998329,     0.22360680},
  {  0.00135425,     0.04034776,     0.36819527,    -0.19277231,     0.13343871, 
     0.04808134,     0.30054689,    -0.29942384,    -0.12281492,    -0.08253271, 
     0.12069781,     0.15301996,     0.55114514,     0.02455466,     0.20698167, 
     0.06560117,    -0.08840462,    -0.04141552,     0.20875959,     0.22360680},
  {  0.03297150,    -0.04486646,     0.01393639,    -0.01230910,     0.00756319, 
     0.00843880,    -0.05882539,    -0.00205577,    -0.03336123,    -0.09705754, 
     0.18223650,     0.27505204,     0.11640776,    -0.98415640,    -0.27634475, 
    -0.25095964,     0.05341233,     0.01884803,    -0.07321042,     0.22360680},
  { -0.00630195,     0.00210376,     0.31845497,    -0.63255608,    -0.47523942, 
    -0.02043996,    -0.22525844,     0.44730283,     0.02518716,     0.03381098, 
     0.05793248,     0.05359087,    -0.12766504,    -0.00338507,     0.09557222, 
     0.16364628,    -0.07974760,    -0.02526895,     0.12321246,     0.22360680},
  {  0.00319372,    -0.02225077,    -0.33584042,     0.36449862,    -0.09373707, 
    -0.00319730,    -0.04373080,     0.29799987,     0.24276207,    -0.25218090, 
     0.10516001,     0.15138497,     0.30068992,     0.01466187,     0.17306533, 
     0.12721931,    -0.07523194,    -0.03266918,     0.17970079,     0.22360680},
  { -0.00169530,    -0.03804556,     0.04368198,     0.03184141,    -0.01882572, 
     0.00999479,    -0.04490470,     0.01325400,    -0.02915776,    -0.07444234, 
     0.13492453,     0.10761661,    -0.14988766,     0.09117883,    -0.03768287, 
    -0.72295676,     0.18139039,    -0.05721761,     0.20834117,     0.22360680},
  { -0.00175328,     0.03765097,     0.06853987,     0.34068518,     0.46048440, 
     0.02471889,     0.29886070,     0.65026178,    -0.76853425,     0.46096107, 
     0.24429166,    -0.15782405,    -0.13606095,    -0.01259603,     0.04012018, 
     0.06745687,    -0.22275879,    -0.03517421,     0.05100764,     0.22360680},
  {  0.67951235,     0.23757404,     0.03081700,     0.05685831,    -0.10137436, 
    -0.02116327,    -0.09565994,    -0.02941115,    -0.19483090,    -0.34342007, 
     0.06052759,    -0.17199625,     0.02087220,     0.02034148,     0.08866969, 
     0.05166274,     0.17972521,     0.12496335,    -0.21266836,     0.22360680},
  { -0.04326895,    -0.03934627,     0.01489937,     0.02459828,     0.01768188, 
    -0.14009872,     0.04628927,    -0.03700050,     0.14049080,     0.35545284, 
     0.22946812,     0.24803819,    -0.01509608,     0.04118764,     0.06436176, 
     0.04688884,     0.11517225,     0.10698573,    -0.24541087,     0.22360680},
  { -0.00768692,     0.06582201,     0.19649492,     0.03710956,     0.49546576, 
    -0.05257424,    -0.47473349,    -0.11358069,     0.10332031,    -0.11619816, 
    -0.03516761,     0.06802488,    -0.33987441,    -0.01860213,     0.15298696, 
     0.19254138,    -0.08722445,    -0.03064029,     0.14288647,     0.22360680},
  { -0.10615431,     0.06920340,    -0.21222304,    -0.14389735,     0.10653252, 
     0.97495615,    -0.14388406,    -0.01288414,     0.08106992,     0.18814040, 
     0.12206715,     0.04417605,    -0.03653889,     0.02169507,     0.08385721, 
     0.05690021,     0.11829105,     0.08840703,    -0.17133399,     0.22360680},
  {  0.00156767,    -0.01618394,     0.00220388,    -0.00511553,     0.00742412, 
     0.02288548,     0.00337393,     0.06659476,    -0.12644854,    -0.16176235, 
    -0.60544385,     0.48698499,     0.04056052,     0.06407434,    -0.08249555, 
    -0.21920208,    -0.35907182,     0.01480767,    -0.37348571,     0.22360680},
  {  0.00065655,    -0.04541786,    -0.01540486,     0.00346833,     0.00192623, 
     0.01169208,    -0.03424510,    -0.05414501,    -0.06613074,    -0.11621025, 
     0.14434250,     0.14036759,     0.02066323,     0.10603287,    -0.84852734, 
     0.30202467,     0.03376821,    -0.00217106,     0.13470635,     0.22360680},
  { -0.17555074,     0.71500757,     0.05252352,     0.15044986,    -0.15168335, 
    -0.01600701,     0.01073625,    -0.04817632,     0.07149679,     0.24853752, 
    -0.23540507,    -0.16128451,     0.03992845,    -0.02678959,    -0.07180661, 
    -0.01917081,     0.02283825,    -0.00444558,     0.09659594,     0.22360680},
  {  0.09808898,    -0.57266478,     0.40288866,     0.36925682,    -0.20178014, 
     0.07748828,    -0.03925209,    -0.03286855,     0.08710476,     0.19608233, 
    -0.26950446,    -0.24739918,     0.03616825,    -0.02368873,    -0.02527910, 
     0.04131159,     0.04546038,     0.02113362,     0.04774680,     0.22360680},
  {  0.00269219,    -0.00277947,     0.00018814,     0.00467786,     0.01060359, 
    -0.00559608,    -0.01301170,     0.00832663,    -0.02661757,    -0.00185008, 
    -0.02040641,    -0.02738710,     0.07616306,     0.00342436,     0.04653034, 
     0.24667282,     0.48175136,    -0.96715232,    -0.61028447,     0.22360680},
  {  0.00655231,    -0.00227730,    -0.00882594,    -0.03739907,    -0.04893686, 
    -0.00941621,    -0.03827905,    -0.14200681,     0.23105441,    -0.04592655, 
     0.45375461,    -0.52670282,     0.01061569,     0.02518422,    -0.13704799, 
    -0.26562166,    -0.64464534,    -0.07863417,    -0.30595534,     0.22360680},
  { -0.66462641,    -0.10786717,    -0.02363327,    -0.05136889,     0.02472298, 
    -0.06458859,     0.01611836,     0.01901677,    -0.13126968,    -0.36582895, 
    -0.02393019,    -0.23873718,     0.02716091,     0.00048965,     0.06083925, 
     0.03374421,     0.16582642,     0.10480953,    -0.15653800,     0.22360680}
 }, {
  { -0.12225778,     0.00643577,    -0.18850007,     0.01297087,    -0.01087376, 
    -0.42269315,    -0.39342375,    -0.05616668,     0.19856921,     0.05036272, 
     0.45650063,     0.05375099,    -0.05381237,    -0.00866985,    -0.01970157, 
     0.07203370,     0.09601553,    -0.04348258,     0.01762634,    -0.22360680},
  {  0.00162841,    -0.00542452,     0.00210746,     0.06206442,    -0.55398519, 
     0.07005425,    -0.26030468,     0.13558391,    -0.20434327,    -0.00473463, 
    -0.29719536,     0.47963384,    -0.08934841,     0.09786738,     0.09787378, 
     0.01274778,    -0.23654047,    -0.10137743,     0.07951741,    -0.22360680},
  { -0.01286989,     0.00304041,    -0.08794279,     0.06754602,    -0.12495226, 
    -0.66125104,     0.60920981,    -0.14756526,    -0.57042931,     0.11041053, 
     0.08271861,    -0.10701010,     0.05091988,     0.02843725,     0.08689631, 
    -0.02038610,    -0.11963992,    -0.15001938,     0.07957345,    -0.22360680},
  { -0.00192516,    -0.00623795,    -0.03264534,    -0.02741543,     0.04372061, 
     0.12140925,    -0.42220748,     0.11025408,    -0.20060114,    -0.03022847, 
    -0.13342195,    -0.59067297,     0.28344793,     0.24391176,     0.11159188, 
     0.00754414,    -0.24552207,    -0.22487113,     0.11566970,    -0.22360680},
  { -0.00700421,    -0.00047631,    -0.00265927,    -0.00027385,    -0.00094533, 
     0.00410918,     0.01325955,    -0.00492586,     0.00478571,    -0.01306550, 
    -0.08041013,    -0.00292185,     0.00211124,     0.02082983,    -0.01782820, 
    -0.02363791,    -0.06522278,    -0.44540653,    -0.95434010,    -0.22360680},
  { -0.08488781,     0.13088507,     0.00703701,    -0.93739912,     0.05651157, 
    -0.01632027,     0.15916606,    -0.09056834,     0.17387555,    -0.00297216, 
    -0.07079224,     0.11574414,    -0.01280400,     0.06399163,     0.06249185, 
     0.01152604,    -0.16459424,    -0.08872350,     0.07009206,    -0.22360680},
  {  0.00763602,    -0.02264113,     0.01055699,     0.22109385,    -0.17578853, 
    -0.03706287,     0.35754272,    -0.03966442,     0.68784332,    -0.08743351, 
    -0.11832652,    -0.06903521,     0.07453004,     0.13635604,     0.05428745, 
     0.03581793,    -0.24253039,    -0.13942052,     0.09534613,    -0.22360680},
  {  0.00337453,    -0.00136488,    -0.01568550,     0.00627964,     0.00271549, 
     0.06946134,     0.03495446,     0.03319039,     0.01050744,    -0.04742082, 
    -0.18865443,     0.00373697,    -0.00821709,    -0.03783935,     0.00931516, 
    -0.10992710,     0.68870836,    -0.32422747,     0.11592903,    -0.22360680},
  {  0.00184283,    -0.00326583,     0.01173330,     0.02418081,     0.01248642, 
     0.03275947,    -0.10328904,     0.06693054,     0.07756103,     0.05985753, 
    -0.11549405,    -0.27034661,    -0.26621118,    -0.70733464,     0.34566297, 
    -0.29114450,    -0.39638555,    -0.10686455,     0.06712878,    -0.22360680},
  { -0.65112644,    -0.11631726,     0.22010293,     0.00343131,     0.00371729, 
     0.01139501,     0.09145038,     0.33422807,    -0.05782687,    -0.21753747, 
    -0.05600829,    -0.09132290,    -0.06040694,    -0.01697076,     0.07377895, 
     0.16823436,     0.10605946,     0.28537202,    -0.05673271,    -0.22360680},
  {  0.02205209,     0.14759051,    -0.00613454,     0.05561145,    -0.00310964, 
     0.02534595,    -0.06590399,    -0.49952044,    -0.00610799,     0.10960731, 
    -0.27662168,    -0.06759333,    -0.01961107,    -0.00739679,     0.05995464, 
     0.10881432,     0.10854761,     0.30126711,    -0.05538419,    -0.22360680},
  {  0.01010280,     0.00070433,     0.00666359,     0.14953082,     0.80010241, 
    -0.04189469,    -0.02962444,     0.06914595,    -0.05322355,    -0.00598735, 
    -0.20165577,     0.34151242,    -0.05146208,     0.09132793,     0.07197956, 
     0.02581761,    -0.21873305,    -0.10947156,     0.07927434,    -0.22360680},
  {  0.47025665,    -0.96282286,    -0.10215314,    -0.18151575,     0.02273412, 
    -0.01147552,    -0.02737209,    -0.26138212,     0.00770170,     0.05768545, 
    -0.10218850,    -0.00933422,    -0.02054646,    -0.00443159,     0.05332822, 
     0.09062080,     0.06557641,     0.20986481,    -0.03379600,    -0.22360680},
  { -0.00377608,     0.01257999,    -0.01419684,    -0.00730488,     0.00024209, 
    -0.02997508,     0.03324683,     0.44135790,     0.10230791,     0.59612537, 
    -0.05095275,     0.12047315,     0.35432975,    -0.05226253,    -0.06290810, 
    -0.23522902,     0.07252686,     0.32605009,    -0.06568819,    -0.22360680},
  { -0.00753582,    -0.00461850,    -0.01142091,    -0.00473590,    -0.01231264, 
     0.00153364,    -0.01076221,     0.01743582,    -0.05219959,    -0.02509269, 
    -0.21422663,    -0.06644671,     0.01043557,    -0.23704162,    -0.80917501, 
     0.16873677,    -0.19094207,    -0.11389295,     0.08423930,    -0.22360680},
  {  0.10730715,     0.00440337,     0.71049525,     0.03221728,    -0.01370892, 
     0.22598514,     0.03710725,    -0.13659240,    -0.06480094,     0.14843475, 
     0.44487094,     0.02593704,    -0.02225666,     0.01590411,    -0.01326667, 
     0.03172478,    -0.02153139,    -0.09476540,     0.04093448,    -0.22360680},
  { -0.13737214,     0.00298731,    -0.62283496,     0.04717431,    -0.00011332, 
     0.54583183,     0.20412460,    -0.06136879,    -0.11660531,     0.10481913, 
     0.41001106,     0.02339159,    -0.04979895,     0.01244994,     0.01082751, 
     0.05586903,    -0.05733325,    -0.02131668,     0.02551377,    -0.22360680},
  { -0.01034308,    -0.00882727,     0.00713715,     0.00855400,     0.00995617, 
    -0.05089080,    -0.00517958,     0.08592925,     0.01168877,     0.07352559, 
    -0.03396331,    -0.34221646,    -0.74928932,     0.57483271,    -0.39702997, 
    -0.74557888,    -0.07605965,     0.33700016,    -0.02708675,    -0.22360680},
  { -0.00064668,     0.00411278,    -0.00088179,     0.01710971,     0.00190976, 
     0.02852769,    -0.02917422,    -0.30468999,    -0.07536053,    -0.67945551, 
     0.22266252,     0.20819045,     0.35654381,    -0.08075816,    -0.08343332, 
    -0.42156182,    -0.07223663,     0.23901977,    -0.02820506,    -0.22360680},
  {  0.54890186,     0.13993203,    -0.05512438,    -0.01802687,    -0.00719042, 
    -0.03015725,     0.05788176,     0.42080197,    -0.02442750,    -0.22781684, 
     0.09069652,    -0.07220345,    -0.07232329,    -0.01354873,     0.05116012, 
     0.15746632,     0.08354236,     0.21184338,    -0.04427932,    -0.22360680},
 }
};

static double d[2][20] = {
 {
  -1.705548e-02,  -1.613788e-02,  -1.437795e-02,  -1.398343e-02,  -1.371158e-02, 
  -1.332744e-02,  -1.244415e-02,  -1.145956e-02,  -1.114267e-02,  -1.043333e-02, 
  -9.509445e-03,  -8.985802e-03,  -7.945402e-03,  -7.161653e-03,  -6.888311e-03, 
  -6.363155e-03,  -6.110669e-03,  -4.959796e-03,  -4.078284e-03,   0.0
 }, {
  -1.815964e-02,  -1.740956e-02,  -1.629862e-02,  -1.544577e-02,  -1.403294e-02, 
  -1.298010e-02,  -1.222090e-02,  -1.212547e-02,  -1.160962e-02,  -1.126248e-02, 
  -9.974435e-03,  -9.118885e-03,  -8.722329e-03,  -8.166752e-03,  -7.612619e-03, 
  -7.089499e-03,  -6.770557e-03,  -5.398941e-03,  -5.000896e-03,   0.0
 }
};

static double f[2][20] = {
{0.077074, 0.050079, 0.046245, 0.053815, 0.014443, 0.040894, 0.063366,
 0.065595, 0.021883, 0.059193, 0.097631, 0.059210, 0.022068, 0.041321,
 0.047703, 0.070746, 0.056779, 0.012674, 0.032359, 0.066921
}, {
 0.07553863, 0.05376433, 0.03768495, 0.04470362, 0.02850415,
 0.03390814, 0.05345819, 0.07803147, 0.03004497, 0.05987184,
 0.09578333, 0.05198703, 0.02191100, 0.04501959, 0.04203491,
 0.06819968, 0.05640919, 0.01573577, 0.03596429, 0.07144495}
};

	for (int i = 0; i < 20; ++i)
	    for (int j = 0; j < 20; ++j)
		a[i][j] = (i == j)? exp(d[m][i] * dt): 0;
	matmpy(a, v[m], a);
	matmpy(a, a, u[m]);
	mattr(a);
#if DEBUG == 2
	printmat(stdout, "%9.6f ", a);
#endif
	vcopy(comp, f[m], 20);
}

/*   NOMALIZATION	*/

static double normlz(double h[][AAS], double fact)
{
	double	trace = 0.;

	for (int i = 0; i < 20; i++) {
	    int	ii = shift_aa(i);
	    for (int j = 0; j <= i; j++)
		h[ii][shift_aa(j)] *= fact;
	    trace += h[ii][ii];
	}
	return (trace / 20);
}

static void makmdm(FILE* fs, FILE* fd, double a[20][20], double comp[20])
{
	double	b[20][20];
	double	c[AAS][AAS];
	double*	nrmf = new double[MAXPAM / PAMSTEP + 1];
	double*	trace = new double[MAXPAM / PAMSTEP + 1];

	for (int i = 0; i < 20; i++) {
	    int	ii = shift_aa(i);
	    for (int j = 0; j < 20; j++)
		c[ii][shift_aa(j)] = b[i][j] = 0.;
	    c[ii][ii] = b[i][i] = 1.;
	}
	int	i = 0;
	for (int pam = 0; pam <= MAXPAM; pam += PAMSTEP, ++i) {
	    double	av, sd;
	    if (pam) logodd(c, b, comp);
	    matstat(&av, &sd, c, comp);
	    nrmf[i] = STDSD / sd;
	    av *= nrmf[i];
	    trace[i] = normlz(c, nrmf[i]);
	    makes(c);
	    printf("Pam = %d  Norl fact = %7.3lf trace = %7.3lf\n",
		pam, nrmf[i], trace[i]);
	    if (fs) putfmtx(fs, c);
	    if (fd) putdmtx(fd, c);
	    matmpy(b, b, a);
	}
	if (fs) {
	    fwrite(nrmf, sizeof(double), i, fs);
	    fwrite(trace, sizeof(double), i, fs);
	}
	if (fd) {
	    fwrite(nrmf, sizeof(double), i, fd);
	    fwrite(trace, sizeof(double), i, fd);
	}
	delete[] nrmf; delete [] trace;
}

static void mkmdmpwr(FILE* fd, double a[20][20])
{
	putpmtx(fd, a);
	for (int i = 1; i < MAXPWR; i++) {
	    matmpy(a, a, a);
	    putpmtx(fd, a);
	}
}

int main(int argc, const char** argv)
{
	int	vt = 0;
	double	a[20][20];
	double	b[20][20];
	double	comp[20];
	double	psetp = (double) PAMSTEP;
const	char*	fn = 0;
static	struct wmode wmd = {1, 1, 0, 0};

	while (--argc > 0) {
	    if (**++argv == '-') {
		switch (argv[0][1]) {
		    case 'c': wmd.tab = argv[0][2] != '-'; break;
		    case 'm': wmd.tab = argv[0][2] != '-'; break;
		    case 'l': wmd.lud = argv[0][2] != '-'; break;
		    case 'p': wmd.pwr = argv[0][2] != '-'; break;
		    case 'f': 
			year = strcmp(argv[0], "-f91")? 1: 0;
			break;
		    case 'v': vt = 1; break;
		    case 'b': vt = 2; break;
		    default:  usage();
		}
	    }
	    else if (**argv == '?') usage();
	    else    fn = *argv;
	}

	if (vt) exp_qdt(a, comp, psetp, vt - 1);	/* vt matrix */
	else {
	    pam1(a, comp);
	    matpwr(a, a, PAMSTEP);	/* a ^ PAMSTEP */
	}
#if DEBUG == 1
	printmat(stdout, "%8.5f ", a);
#endif
	matcpy(b, a);
	for (int i = 0; i < 20; i++) {
	    if (i % 10 == 0) putchar('\n');
	    printf(" %7.4f", 100. * comp[i]);
	}
	putchar('\n');

	FILE*	fs = 0;
	FILE*	fd = 0;
	if (wmd.cmp) {
	    fd = fopenpbe(fn, mdm_cmp, 0, "wb", 1);
	    fwrite(comp, sizeof(double), 20, fd);
	    fclose(fd);
	    fd = 0;
	}
	if (wmd.tab) fs = fopenpbe(fn, mdm_tab, 0, "wb", 1);
	if (wmd.lud) fd = fopenpbe(fn, mdmld_tab, 0, "wb", 1);
	makmdm(fs, fd, a, comp);
	if (fs) fclose(fs);
	if (fd) fclose(fd);
	if (wmd.pwr) {
	    fd = fopenpbe(fn, mdmpwr_tab, 0, "wb", 1);
	    mkmdmpwr(fd, b);
	    fclose(fd);
	}
	return(0);
}
